library(spdep)

random_submap <- function(geo_obj, n_units){

    al <- poly2nb(geo_obj, queen = FALSE)
    
    ## Sample seed
    ds <- sample(1:length(al), 1)
    visited <- c()

    while(length(ds) < n_units){
        
        ## Sample node to check neighbors of
        get_adj <- as.numeric(sample(as.character(ds[!(ds %in% visited)]), 1))
        candidates <- sample(al[[get_adj]])
        candidates <- candidates[!(candidates %in% ds)]
        
        ## Loop through neighbors, add
        if(length(candidates) > 0){
            for(i in 1:length(candidates)){
                ds <- c(ds, candidates[i])
                visited <- c(visited, get_adj)
                if(length(ds) == n_units){
                    break
                }
            }
        }
        
    }

    ## Return map
    return(geo_obj[ds,])

}
